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1 Introduction 


An important area of modern research is global optimization as it occurs very frequently 
in applications (extensive surveys on global optimization can be found in Neumaier 
fl^ . Floudas [3,111, Hansen [i1|i and Kearfott [ill). A method for solving global 
optimization problems efficiently is by using a branch and bound algorithm (^ e.g ., 
BARON by Sahinidis I 2 I. S . the COCONUT environment by Schichl [jj. lil. [ia . 
or LINGO by Schrage l27ll L which divides the feasible set into smaller regions and 
then tries to exclude regions that cannot contain a global optimizer. Therefore, it is 
important to have tools which allow to identify such regions. In this paper we will 
present a method which is able to find such regions for a CSP (constraint satisfaction 
problem), i.e. for a global optimization moblem with a constant objective function, by 
generalizing the approach from Fendl [^. 

Certificate of infeasibility. For this purpose we consider the CSP 


F { x ) G F 

X € X 


( 1 ) 
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with F : K” —>■ M™, x € IM", F € IM™, and we assume that a solver, which is able 
to solve a CSP, takes the box u := [u,u] C x into consideration during the solution 
process. We constructed a certificate of infeasibility /, which is a nondifferentiable and 
nonconvex function in general, with the following property: If there exists a vector y 
with 

f{y,u,u)<0, (2) 

then the CSP IH) has no feasible point in u and consequently this box can be excluded 
for the rest of the solution process. Therefore, a box u for which 0 holds is called an 

exclusion box. 

Easy examples immediately show that there exist CSPs which have boxes that 
satisfy ( 0 , so it is worth to pursue this approach further. 

Exclusion boxes. The obvious way for finding an exclusion box for the CSP o is to 
minimize / 


mm/{y,u,'u) ( 3 ) 

and stop the minimization if a negative function value occurs. Since modern solvers 
offer many other possibilities for treating a box, we do not want to spend too much 
time for this minimization problem. Therefore, the Idea is to let a nonsmooth solver 
only perform a few steps for solving ©• 

To find at least an exclusion box v := [u, tJ] C u with y + r < v, where r £ (0, u —u) 
is fixed, we can try to solve the linearly constrained problem 

min f{y,y, v) 

y,2l,v 

s.t. [u + r, tJ] C li . 

Another important aspect in this context is to enlarge an exclusion box v by solving 

max/r(w, v) 

V,v,v (-4) 

S.t. f{y,y,v) < S , \y,v] C u , 

where <5 < 0 is given and y, measures the magnitude of the box v (e.g., y{y_,v) : = 
1^ “ nil)- Since only feasible points of are useful for enlarging an exclusion box 
and we only want to perform a few steps of a nonsmooth solver as before, we expect 
benefits from a nonsmooth solver that only creates feasible iterates because then the 
current best point can alwws be used for our purpose. For proofs in explicit detail we 
refer the reader to Fendl [J p. 147 ff. Chapter 5]. 

The paper is organized as follows: In Section [2] we first recall the basic facts of 
interval analysis which are necessary for introducing the certificate of infeasibility which 
is done afterwards. Then we discuss some important properties of the certificate and 
we explain in detail how the certificate is used for obtaining exclusion boxes in a CSP 
by applying a nonsmooth solver. In Section[3]we explain how we obtain a starting point 
for optimization problems to which we apply the nonsmooth solver. 

Throughout the paper we use the following notation: We denote the non-negative 
real numbers by ]R>o := {x € M : a: > 0 } (and analogously for < as well as >). 
Furthermore, we denote the p-norm of a; £ K” by \x\^ for p £ {1, 2, 00}. 
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2 Presentation of the application 

After summarizing the most brisic facts of interval analysis, we construct the certificate 
of infeasibility in this section. Furthermore, we discuss how a nonsmooth solver can 
use this certificate to obtain an exclusion box in a CSP. 


2.1 Interval arithmetic 

We recall some basic facts on interval arithmetic from, e.g., Neumaier [Till . We denote 
a box (also called interval vector) by a: = [ x , x ] and the set of all boxes by IK” := {x : 
X = [ x , x \, X < x } . For S' C R bounded, we denote the hull of S by DS := [inf S, sup S]. 
We extend the arithmetic operations and functions ip : R —> R to boxes by defining 

p>{x) := □{(p(a;) : x € a;} . (5) 

For every expression of (p : R” —y R which is a composition of arithmetic operations 
and elementary functions the fundamental theorem of interval arithmetic holds 

[inf <p(x),sup(p(x)] C ${x) . (6) 

2:6a: xGa; 


2.2 Certificate of infeasibility 

Let s : R"* x R” x R'^ x IR” —v IR be a function. We assume that Z : R'” x R” x 
R« X R” X R” —^ R 

Z{y,z,w,x,x) :=s\i^s{y,z,'w,x) , (7) 

where x = [x, x] G IR”, satisfies 

Z{y, z, w, X, x) > sup (^{x) - F(z)) . (8) 

xGx 

Example 1 If we set w = {R, S) G R”j.fJ^ x R”^^;”, where we denote the linear space of 
the upper resp. strictly upper triangular n x n-matrices by R”j.fJ* — R”^ with ni : = 
^n(n + 1 ) resp. R”tji” — R”“ with no := ^{n — l)n, which implies q = no + ni = , 

and if we define 

m 

si{y,z,R,S,x) := (^'^ykdk[z,x] + {x - z)'^{R^R + s'^ - S)^{x - z) (9) 

k=l 

then the corresponding Z satisfies ([51) because: Due to the skew-symmetry of — S, 
we have 

y'^{F{x)-F{z))+{x-zf{R^R + s'^-S){x-z)>y'^{F{x)-F{z)) . (10) 

for all X G a:. Since the slope expansion 


Fk{x) = Fk{z) + Fk[z,x]{x - z) , 


(11) 
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where the slope Fi;[z,x] € holds for all x,z £ M” (cf., e.g., Neumaier El) , we 

obtain for all rr E aj 

y^{F{x) - F{z)) F{x- zf{R^R + - S){x - z) C si{y, z, R, S, x) (12) 

due to ([TT1). ©, (O, and ®- Now we obtain ((81) due to ([Toll, (dll), and 0. 


Proposition 1 It holds for all z £ x 

Z{y, z,w,x,x) > 0 . (13) 

Proof. (I13II follows from (|H|) and the assumption that z £ x. □ 

M and / : K™ X K" X ]R« X K" X P” 


Definition 1 We define Y : R"* x M" 
by 


Y{y,z) ■=mIy'^{F-F{z)) 

Z{y,z,w,x,x)-inaic(0,Y{y,z)) 

f{y,z,w,x,x) := - — -^ 

T[y,z,w,x,x) 


(14) 

(15) 


where T : R™ x K" x R® x R" x R’" 
almost everywhere. 


R>0 is positive, continuous and differentiable 


Remark 1 f from m depends on N = m + 3n + q variables and is not differentiable 
everywhere and not convex (in general). 


Now we state the main theorem for our application. 

Theorem 1 If there exist y £ R™, x<z<x£ R", and w £M.‘^ with f{y, z, w, x, x) < 
0, then for all x £ x there exists k £ {1,.. ., m} with F/^(x) 0 F/^, i.e. there is no x £ x 
with F{x) £ F, i.e. there is no feasible point. 

Proof, (by contradiction) Suppose that there exists x £ x := [x,x] with F{x) £ F. 
By assumption there exist y £ R"* and z £ x C IR” with f{y,z,w,x,x) < 0, which is 
equivalent to Z{y,z,w,x,x) < Y{y,z) due to (11511 and (1131) . Since 

Z{y,z,w,x,x) > snpy'^[F{x) - F{z)) > y'^{F{x) - F{z)) 

X^X 


for all a; € a: due to ((8| and 

Y{y,z)< inf y^{F-F{z)) <y^{F-Fiz)) 

F€F 

for all due to (jldp and ([21), we obtain for all x ^ a? and for all 

y'^ [F{x) - F{z)) < Z{y,z,w,x,x) < Y{y,z) < y'^ [F - F{z)) , 

which implies that we have y'^F{x) < y^F for all a; G a: and for all F £ F. Now, 
choosing x = x £ x and F = F(x) G P in the last inequality yields a contradiction. □ 

The following proposition gives in particular a hint how the y-component of a 
starting point should be chosen (cf. (1331) 1. 
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Proposition 2 We have for all y G and z G M" 


Y{y,z) 


f VkiEk - Fkiz)) for 1/fe > 0 
^ 1 yk(Fk - Fk{z)) for yk<0 . 


(16) 


Furthermore, let I, J C {1,..., m} satisfy / V J ^ F_^ = —oo A yi > 0 for all 

i € I and F j = oo A yj < 0 for all j G J, then we have for all z G K" 

Y{y,z) = -oo. (17) 


Proof. 

F, = 

(HSJ. 


US holds because of (O- For obtaining (nzi, consider without loss of generality 
—oo A > 0 and F2 = oo A 2/2 < 0, then the desired result follows from 

□ 


2.3 Properties of the certificate for quadratic F 

In this subsection we consider the special case of / with quadratic F, i.e. 

Fk{x) = clx + CkX (18) 

with Cfc G K” and Cfc G for k = . ,m. Since 

Fk{x) - Fk{z) = {cl + x^Ck + z^cl) {x - z) 
for all x,z £ M", the slope expansion from m holds with 

Fk[z,x] = cl + x^Ck + z'^cl . (19) 


Proposition 3 Let Fk be quadratic and set 

m m 

C{y) ■='^Ckyk, c(y,z) ■.= '^Ckyk + {C{y) + C{y)'^)z, A{y,R,S) 

( 20 ) 

■.= C{y) + R^R + s'^ -S . 

Then ((HI) is satisfied for 

S 2 {y,z,R,S,x) := {c{y,z)'^ + {x - z)'^A{y, R, S)){x - z) . (21) 

Proof. Since HkliykFklzjX] = c{y,z)'^ + {x — z)^C{y) due to (fT^ and (l20ll . we 
obtain Ylk=i ykFk[z,x] + (a; - z)'^{R^R + S'^ - S) = c{y,z)'^ + (a; - z)'^A{y,R,S) 
due to (|20l). which implies that S 2 from mi) has the same structure as si from 
and consequently we obtain that ((SJ holds for S 2 , too. □ 


Proposition 4 Let Fk be quadratic, then we have for all p G M"* and q G K 

C{y + ap) = C{y) + aC{p) , c{y + ap, z) = c{y, z) + ac{p, z) (22) 


and furthermore we have for all k, > 0 

A{K^y, kR, S) = A{y, R, S) . 


(23) 
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Proof. (I22| holds due to dSOj. 1231) holds due to 1201) and HH). □ 

Proposition 5 Let be quadratic. If the positive function T : M"* x K" x x 
K”” X M" X R" —^ R>o satisfies the (partial) homogeneity condition 

T{K^y, z, kR, K^S,x,x) = K'^T{y,z,R,S,x,x) (24) 

for all K, > 0, y £ R"*, 2 £ [x,x] £ IR", R £ and S £ R"triu> ^^6 certificate 

f from US) is (partially) homogeneous 


f[r?y,z,nR,r?S,x,x) = f{y,z,R,S,x,x) . (25) 

Proof. Since F}. is quadratic by assumption, the statements of Proposition [3] hold. By 
using 1221 ) and 123) we calculate 

c ( k ^ i /, 2 )^ + (a 3 - z)^A{hi^y, kR, r?S) = [c{y, 2 )^ + {x - z)^A{y, R, S)) (26) 

and therefore Z{K^y,z,KR,K^S,x,x) = Z[y,z,R,S,x,x) follows due to Q, (1^ . 
and 126|. Furthermore, we have Y{n?'y,z) = K^Y{y,z) due to (1141) and, hence, we 
obtain max (0,Y(K?y, z)') = max (O, F(i/, 2 )). Consequently, (I15II and (I24II imply 

(E51) . □ 

Remark 2 The intention of 125]) is to reduce the scale dependence of the unbounded 
variables y, R and S of /. If we go through the proof of Proposition IS] again, we notice 
that we use the scaling property (1231) of A for showing (1261) . From the proof of (1231) we 
notice that this proof only holds, if y, R and S are treated as variables and none of 
them is treated as a constant (since factoring out of a constant, yields an additional 
factor K~^ to the constant). Nevertheless, if one of the variables y, R resp. S is treated 
as a constant and we set the corresponding value to 1 / = 0, i? = 0 resp. S = 0, then 
the proof still holds. 

Example 2 Consider the variables R and S as constants and set R = S = Q. Then 
Ti{y,z,R,S,x,x) := 1 does not satisfy (1^ . while T 2 {y,z,R,S,x,x) := I 2 /I 2 does 
(cf. Remark 121). Note that T 2 violates the requirement of positivity, as demanded in 
Definition ITJ for y = 0, and hence in this case / is only defined outside the zero set 
of T. Nevertheless, since the zero set of T 2 has measure zero, it is numerically very 
unlikely to end up at a point of this zero set and therefore we will also consider this 
choice of T due to the important fact of the reduction of the scale dependence of / as 
mentioned in Remark 12] (also cf. Subsection 12.41 (directly after optimization problem 
(1271) 1 and Example 13 . 

Example 3 Choose n = 2, m = 2, ci = ( _^3 ) , C 2 = ( 2 ), C*! = ( 3 4 ): ^2 = ( Z 2 7 ) ’ 
which yields Fi{x) = 2x\ +xi + 3 a;ia :2 — 3*2 + 4*2 and + 2 ( 2 :) = —*1 +4*i — 2 * 1*2 + 
2*2 + 7*2 due to (1181) . x = [—3,3] x [—4,4] and F = [—1,7] x [—2,0], then we can 
illustrate certificate / from 113 with different T from Example [2] by the following plots 
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2.4 Exclusion boxes for constraint satisfaction problems 

Now we explain in detail how to use Theorem[T]for finding exclusion boxes for the CSP 
©• If we apply a solver for linearly constrained nonsmooth optimization (Note: The 
certificate / from m is not differentiable everywhere due to Remark [T]) to 

min/{i/, 2 , w,M,n) 

s.t. u,v £x , u + r <v , z £[u,v\ (27) 

y € K"*, z,u,v € R", U) G 

with a fixed r £ [0,x — — although the convergence theory of many solvers (cf., 

e.g., the second order bundle algorithm by Fendl & Schichl [3, p. 7, 3.1 Theoretical 
basics]) requires that all occurring functions are defined on the whole R^, which might 
be violated for certain choices of T (cf. Example [2]) — and if there occurs a function 
value smaller than zero (during the optimization process), then there is no feasible 
point in [m, v] according to Theorem [T] and consequently we can reduce the box x to 
the set X \ [u, n] in the CSP ([T|). 

If [u, v] = X (i.e. u and v are fixed and therefore no variables), then we can reduce 
the box X to the empty set, i.e. the reduction of a box to the empty set is equivalent 
to removing the box. 

The constant r determines the size of the box [u, v], which should be excluded: The 
closer r is to 0, the smaller the box [m, u] can become (if r = 0, [u, v] can become thin, 
what we want to prevent, since we want to remove a preferably large box [m, w] out of 
X, as then the remaining set x \ [m, n] is preferably small). 

If during the optimization a point z £ x is found with F{z) G F, then we have 
found a feasible point and therefore we can stop the optimization, since then we cannot 
verify infeasibility for the box x. 

Remark 3 If y C cc and we remove y from x, then the remaining set x\y is not closed. 
Nevertheless, if we just remove y° C y, then the remaining set x\y° D a: \ y is closed 
(i.e. we remove a smaller box and therefore the remaining set is a bit larger, since it 
contains the boundary of y). Furthermore, the set x\y can be represented as a union 
of at most 2n n-dimensional boxes, i.e. in particular the number of boxes obtained by 
this splitting process is linear in n. 

We make the assumption that the certificate of infeasibility from (US of the box 
[m, D] satisfies /(y, z, w, u, v) =: 5 < 0, i.e. [m, D] is an exclusion box according to 
Theorem [T] For <5 G [5, 0) and a box x with [u, D] C x, we can try to apply a solver for 
nonlinearly constrained nonsmooth optimization to 

minti(y, z, w, u, v) 

s.t. f{y, z, w, u.v) < 5 

7 . (28) 
U,V (z X , U < U , V < V , z G [u,v] 

y G R™, z,u,v £ R", u; G R'^ 

to enlarge the exclusion box [rt, w] in x, where b : R^ —> R is a measure for the box 
[u, v\ in the following sense: If b ■ rN -^ 

, then the following conditions must hold: 
If [u, n] is small, then 6(., u, v) is close to 0 and if [m, n] is large, then b(., u, v) is negative 
and large. This means: The larger the box [u, i>] is, the more negative b{.,u,v) must 
be. For examples of this type of box measure cf. (1301) . Alternatively, if b : R^ ^R>o, 
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then the following condition must hold: If [u, u] C a: is close to x, then b{.,u, v) is close 
to 0. For examples of this type of box measure cf. m- 

Remark 4 In opposite to (123), where the linear constraint u + r < v occurs, we use in 
dm the bound constraints u <u and v < v. 

Furthermore, we make the following two observations for the global optimization 
problem 


minFobjla:) 

X 

s.t. Fi{x) € Fi for alH = 2 ,..., m (^9) 

X £ X , 

where Fobj : K R: First of all, the certificate / from m can be used for finding 
exclusion boxes in the global optimization problem 11291) with an arbitrary objective 
function Fopj, since the certificate / only depends on the constraint data F, F and x 
(cf. the CSP ([3) and since a solution of an optimization problem is necessarily feasi¬ 
ble. Secondly, we denote the current lowest known function value of the optimization 
problem (1291) by Now, if we can find a box [u,u] C x for which the certificate / 

from (1151) with Fi := [—oo,FQ((f] has a negative value, then Theorem [T] implies that 
for all X £ [u, u] there exists k £ {1,..., m} with F/^(x) ^ F/^, which is equivalent that 
for all X £ [m, u] we have < ^ 1 ( 3 ;) or there exists i £ {2,... ,m} with Fi{x) ^ Fi, 
i.e. any point in the box [u, w] has an objective function value which is higher than the 
current lowest known function value ^^((1 or is infeasible. Consequently, the box [u, v] 
cannot contain a feasible point with function value lower equal and hence the 

box [u,v\ cannot contain a global minimizer of the global optimization problem (1291) . 
Therefore we can exclude the box [u, u] from further consideration. 

Example 4 For measuring the box [u, u] with b : —>■ R<o, we can use any negative 

p-norm with p £ [1, 00 ] as well as variants of them 

b\{y,z,'w,u,v) := -|w -uli , b‘i{y,z,w,u,v) 

b’^{y,z,w,u,v) := -lu-uL ■ 

For measuring the box [u, u] with b : R^ — R>o, we can use any p-norm with 
p € [I, cx)] as well as variants of them 

b\{y,z,w,u,v) . b\{y,z,w,u,v) 

b+{y,z,w,u,v) := |(^“Z|)L ■ 

b- is concave, while b+ is convex. b_ can be used for unbounded x, while this is not 
possible for is smooth, while b^ and b°° are not differentiable, b^ has an equal 

growing rate for all components. The growing rate of b^ depends on sgn (^|.|2 — l). 
already grows, if the absolute value of the largest components grows. 

Example 5 Choose n = 1, m = I, w = 0, ci = 1, Ci = ^, which yields F{x) = ^x^ + x 
due to (HU), as well as a: = [—1,2] and consider two CSPs ((3 with F^ = [—2,1] 
resp. F^ = [—2, —1] which yield the following two graphics 



1 I 12 

:= -Trk; - uL 


(30) 



10 


Hannes Fendl et al. 




from which we can see that the CSP has feasible points for F^, while it is infeasible 
for F'^. The corresponding certificates / from GSl), where we only consider the vari¬ 
ables y £ K and z G a: as well as different T from Example [5] and where we denote 
the function value of a local minimizer of the optimization problem (^7)) by /, can be 
illustrated by the following plots 




Fig. 13; / for {F^,T 2 ) ^ / = 0 Fig. 14: / for {F^,T 2 ) ^ / = -| . 

We see from Figure [T51 and Figure [T^ that the certificate / is not defined for y = 0 due 
to the definition of T 2 in Example [2] 
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3 Starting point 


We implemented the suggestions from Subsection l2.4l in GloptLab by Domes which 
is a configurable MATLAB framework for computing the global solution of a quadratic 
CSP ([T|), i.e. with Pfc from (11811 . The matrices Ck € are lower triangular in 

GloptLab 

Gfe G i" := {A G M”''” : = 0 for i < j} . (32) 

For running GloptLab the MATLAB toolbox INTLAB by Rump lilll . lp_solve by 
Berkelaar et al. 0, SEDUMI by Polik [^, Sturm 0 as well as SDPT3 by Toh 
et al. (2^ were installed for using all features of GloptLab. 

So the last issue that remains to be discussed is, how to find a point {y, z, w, u, v) be¬ 
ing feasible for the linearly constrained optimization problem (I27II with f{y, z, w, u, v) < 
0 quickly. For this we need a good starting point (y^, z ^and therefore we 
must take the following observations into account: y^ and z^ should be chosen so, 
that Y[y'^,z^) is positive, and [y'^, z'^,w^,u^,v^) should be chosen so, that the term 
Z[y'^, z^,nP,vP,v^), which is non-negative due to (11311 . is near zero. These facts lead 
to the following suggestions for choosing a starting point [y'^, z^ ,w^ ,vP ,v^)\ First of 
all, if the solver in use can only handle strictly feasible bound/linear constraints (e.g., 
the second order bundle algorithm by Fendl & Schichl 0 with using socp by Lobo 
et al. 0 for computing the search direction), then the initial choices of vP and 
must satisfy vP,v^ G {x,x) and + r < v^, e.g., := (1 — to)* + to(* ~ *) and 

:= (1 — ti)^ + ti{x — r) for some fixed toTi ^ (0,1) with to < fi- Otherwise (e.g., 
SolvOpt by Kappel & Kuntsevich fllll or the second order bundle algorithm with 
using MOSEK for computing the search direction) we take the endpoints of x for vP 
and v^. Secondly, the natural choice for the starting value of z G [w, n] C a: is the mid¬ 
point := ^(m*^-|- w*^) of the box Thirdly, to get the term max (O, T(y, z)) in 

the certificate / from m as large as possible, we make the following choices: For the 
case = —00 resp. the case Ej, = 00 resp. the case that both and E^ are finite, 
we choose 


0.^ f-lifEfe<Efe(z°) 0 .^ flifEfc(zO)<Efe 0 
\ 0 else \ 0 else ’ 


lifEfc(z°)<E, 
-lifEfe<Efc(zO) 
0 else 

(33) 


respectively due to m and HZD- Finally, for the choices of R and S we refer to 
Proposition 


Remark 5 If we choose T = T 2 (cf. Example [21, then T{y,z, w,u,v) = 0 <=> y = Q. 
Therefore, if = 0 occurs as starting point, then we have a feasible point E(z) G F 
due to (|33l). Furthermore, we can expect that no solver should have difficulties with 
this choice of T because of the small size of the zero set of T due to Example (21 


In the following we will make use of the MATLAB operators diag, tril and triu. 

Proposition 6 Let Ej, be quadratic and let (1321) be satisfied. Choose any y G and 
consider the modified Cholesky factorization 


A = fFr - D 


( 34 ) 
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of A (with R € R”j.hr non-negative diagonal matrix D £ where 

i := C{y) + - S G , S := -\Xrm{C{yf, l) G K^triu (35) 

and denotes the space of all symmetric n x n-matrices. Then 

A = C(y) - itril(C{i/), -l) + ^tnu(C(y)^, l) , diag(A) = diag(C'(3/)) . (36) 

Furthermore, if we set 

R:=Di , (37) 

then A{y, R, S) = ff’R G . 

Proof. Since is quadratic by assumption, the statements of Proposition[3]hold. Since 
C{y) is lower triangular due to (I32II and (I20II . we obtain S G R”triu 

A = C{y)-^^tTm(C{y)'^,l)'j + itriu(C'(i/)^, l) (38) 

due to (1551) . Now, (13811 implies (1361) . We calculate 

^(y)^ - 5triu(C(y)'^, l) = diag(C(y)) + itriu(C(y)'^, l) . (39) 

Because of ^ (triu(C(i/)^, l))^ = tril(C'(i/), —l) —^ (triu(C(y)'^, l))^ we have diag(C'(j/)) + 
i (triu(C'(j/)^, l))^ = C(2/) — ^ (triu(C'(t/)^, l))^. Therefore, combining ll38ll and (l39ll 
yields A^ = A, i.e. A G R”y^. Consequently there exists a modified Cholesky factor¬ 
ization of A of the form (13411 . Hence, we can choose R according to (13711 and evaluating 
A at {y, R, S) with R from (13711 and S from (13511 yields A{y, R, S) = R^R due to (I2UI1 . 

([37)1 and ((Ml)- □ 

Remark 6 If T is positive semidefinite, then D = 0 due to (13411 . If C{y) is a diagonal 
matrix, then 5 = 0 due to (135)1. Due to dMI, we can construct A by setting A equal 
to C{y), then multiplying the lower triangular part of T by ^ and finally copying the 
resulting lower triangular part of A to the upper triangular part of A. 

Now we combine the facts that we presented in this subsection to the Algorithm [2l 
which we will use for creating a starting point for the optimization problem (12711 with 
quadratic F. 

Algorithm 2. 

if the solver can only handle strictly feasible bound/linear constraints 
Choose 0 < to < ti < 1 (e.g., to = 0.1 and t\ = 0.9j 
u = (1 - to)x -h to(T - r) 
u = (1 - ti){x-fr) 4- tfx 
else 

u = X 
V = x 

end 

Z = ^(u-hw) 

F = F{z) 
if F € F 

stop (found feasible point => cannot verify infeasibility) 
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end if 

for k = 1 ■. m 

£fc = -oo 
if Fk < Fk 
Vk = -1 

else 

yk = o 

end if 

else if = oo 
if Fk < F{k) 

Vk = 1 
else 

yk = o 

end if 
else 

if Fk < Fk 
t/fc = ^ 

else if Fk < Fk 
yk = -1 
else 

yk = o 

end if 
end if 
end for 
Compute C(y) 

S = -ftriu{C{y)'^, l) 

[R,D] =modified_cholesky_factorization(^C{y) + S’^ — S) 
R =sqrt(^D) 


Remark 7 Infeasible constrained solvers (e.g., SolvOpt by Kappel & Kuntsevich 
0) can be applied directly to the nonlinearly constrained optimization problems (1281) . 
In this case the starting point created by Algorithm [5] can be used at once without 
solving optimization problem Il27llfi rst as it is necessary for the second order bundle 
algorithm by Fendl & Schichl Q. Therefore, the bound constraints u < u, v < v oi 
optimization problem (12811 do not occur in this situation. Nevertheless, it is useful in 
this case to add the linear constraint u + r < v (with a fixed r > 0) from optimization 
problem m to the constrained problem for preventing the box [u, u] from becoming 
too small. 


4 Numerical results 

In the following section we compare the numerical results of the second order bundle 
algorithm by Fendl & Schichl [^, MPBNGC by MAkela fl^ and SolvOpt by 
Kappel & Kuntsevich [13 for some examples that arise in the context of finding 
exclusion boxes for a quadratic CSP in GloptLab by Domes [^. 
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4.1 Introduction 


We will make tests for 


(the reduced version of) the second order bundle algorithm for nonsmooth, non- 
convex optimization problems with inequality constraints by Fendl & Schichl 
6 | (with optimality tolerance e := 10~® and with MOSEK by Andersen et al. 

as QCQP-solver for determining the search direction), where we refer to the 
linearly constrained version as “BNLC” and to the nonlinearly constrained version 
as “Red(uced) Alg(orithm)”. It is an extension of the bundle-Newton method for 
nonsmooth, nonconvex unconstrained minimization by Luksan & Vlcek fl^. [i3 
to the nonlinearly constrained problems. 

MPBNGC by Makela fl^ (with standard termination criterions; since MPBNGC 
turned out to be very fast with respect to pure solving time for the low dimensional 
examples in the case of successful termination with a stationary point, the number 
of iterations and function evaluations was chosen in a way that in the other case 
the solving times of the different algorithms have approximately at least the same 
magnitude) 

SolvOpt by Kappel & Kuntsevich [llll (with the standard termination criterions, 
which are described in Kuntsevich & Kappel [i3| I 


(we choose MPBNGC and SolvOpt for our comparisons, since both are written in a 
compiled programming language, both are publicly available, and both support non¬ 
convex constraints) on the following examples: 

— We give results for the linearly constrained optimization problem (p7l) with a fixed 
box (i.e. without optimizing u and v) for dimensions between 4 and II in Subsection 

S31 

— We give results for the linearly constrained optimization problem II27I) with a vari¬ 
able box (i.e. with optimizing u and v) for dimensions between 8 and 21 in Sub¬ 
section SZl 

— We give results for the nonlinearly constrained optimization problem (1281) for di¬ 
mension 8 in Subsection [431 where we use b]_(y, z, R, S, u, v) := \ ^“_=^ as the 
objective function. 


The underlying data for these nonsmooth optimization problems was extracted from 
real CSPs that occur in GloptLab by Domes [^. Apart from u and v, we will concen¬ 
trate on the optimization of the variables y and z due to the large number of tested 
examples (cf. Subsection 14.211 . and since the additional optimization of R and S did 
not have much impact on the quality of the results which was discovered in additional 
empirical observations, where a detailed analysis of these observations goes beyond the 
scope of this paper. Furthermore, we will make our tests for the two different choices 
of the function T from Example (2] which occurs in the denominator of the certificate / 
from (O, where for the latter one / is only defined outside of the zero set of T which 
has measure zero. 

The (extensive) tables corresponding to the results, which we will discuss in this 
section, can be found in Fendl & Schichl Q. 

All test examples will be sorted with respect to the problem dimension (beginning 
with the smallest). Furthermore, we use analytic derivative information for all occurring 
functions (Note: Implementing analytic derivative information for the certificate from 
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(I15II effectively, is a nontrivial task) and we perform all tests on an Intel Pentium IV 
with 3 GHz and 1 GB RAM running Microsoft Windows XP. 

We introduce the following notation for the record of the solution process of an 
algorithm. 

Notation 3. We denote the number of performed iterations by Nit, we denote the final 
number of evaluations of function dependent data by 

Na := “Number of calls to (/, g, G, F, g, G)” (Red Alg) 

Nb := “Number of calls to (/, g, F, g)” (MPBNGC) 

Nc := “Number of calls to (/, F)" (SolvOpt) 

Ng := “Number of calls to g” (SolvOpt) 

Ng := “Number of calls to g” (SolvOpt) , 

and we denote the duration of the solution process by 

ti := “Time in milliseconds” 

t2 := “Time in milliseconds (without (QC)QP)” , 

where t2 is only relevant for the second order bundle algorithm . 

Remark 8 In particular the percentage of the time spent in the (QG)QP in the second 
order bundle algorithm is given by 

_ ti(Red Alg)-t2(Red Alg) / 

Pi ■— ti(Red Alg) ■ 

For comparing the cost of evaluating function dependent data (like e.g., function 
values, subgradients,...) in a preferably fair way (especially for solvers that use different 
function dependent data), we will make use of the following realistic “credit point 
system” that an optimal implementation of algorithmic differentiation in backward 
mode suggests (cf. Griewank & Corliss 0 and Schichl [23.12^.12^1. 

Definition 2 Let /a, gA £did Ga resp. Fa, (Ja ^md Ga be tbe number of function 
values, subgradients and (substitutes of) Hessians of the objective function resp. the 
constraint that an algorithm A used for solving a nonsmooth optimization problem 
which may have linear constraints and at most one single nonsmooth nonlinear con¬ 
straint. Then we define the cost of these evaluations by 

c{A) := fA + 3gA + 3 N-Ga + nlc • {Fa + 3gA + 3N ■ Ga) , (41) 

where nlc = 1 if the optimization problem has a nonsmooth nonlinear constraint, and 
nlc = 0 otherwise. 

Since the the second order bundle algorithm evaluates /, g, G and F, g, G at every 
call that computes function dependent data (cf. Fendl & Schichl Q), we obtain 

c(Red Alg) = (1 + nlc) • Na • (1 -h 3 -h 3N) . 

Since MPBNGC evaluates /, g and F, g at every call that computes function dependent 
data (cf. MakelA 0), the only difference to the second order bundle algorithm with 
respect to c from su is that MPBNGC uses no information of Hessians and hence we 
obtain 

c(MPBNGC) = (1 + nlc) • Nb • (1 -h 3) . 
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Since SolvOpt evaluates / and F at every call that computes function dependent data 
and only sometimes g or g (cf. Kuntsevich & Kappel [l^ l. we obtain 

c(SolvOpt) = (1 + nlc) • Nc + 3{Ng + nlc • Ng) . 

We will visualize the performance of two algorithms A and B in this section by 
the following record-plot: In this plot the abscissa is labeled by the name of the test 
example and the value of the ordinate is given by rp(c) := c{B) — c{A) (i.e. if rp(c) > 0, 
then rp(c) tells us how much better algorithm A is than algorithm B with respect to 
c for the considered example in absolute numbers; if rp{c) < 0, then rp(c) quantifies 
the advantage of algorithm B in comparison to algorithm A-, if rp(c) = 0, then both 
algorithms are equally good with respect to c). The scaling of the plots is chosen in 
a way that plots that contain the same test examples are comparable (although the 
plots may have been generated by results from different algorithms). 


4.2 Overview of the results 


We compare the total time of the solution process 



ti (Red Alg) 

t'p (Red Alg) 

PI 

ti (MPBNGC) 

ti (SolvOpt) 

Linearly constrained (fixed box) 

1477 

215 

T = 1 

0.85 

231 

2754 

Linearly constrained (variable box) 

782 

60 

0.92 

30 

1546 

Nonlinearly constrained 

25420 

4885 

0.81 

21860 

38761 

Nonlinearly constrained (*) 

19053 

3723 

0.80 

2067 

30312 

Linearly constrained (fixed box) 

1316 

129 

T = |y|2 
0.90 

15 

1508 

Linearly constrained (variable box) 

797 

45 

0.94 

30 

2263 

Nonlinearly constrained 

24055 

4284 

0.82 

25383 

16909 

Nonlinearly constrained (*) 

18038 

3112 

0.83 

3719 

12635 


where we make use of (SSI and in (*) we consider only those examples for which 
MPBNGC satisfied one of its termination criterions (cf. Subsection 14.511 . 

For the linearly constrained problems MPBNGC was the fastest of the tested al¬ 
gorithms, followed by BNLC and SolvOpt. If we consider only those nonlinearly con¬ 
strained examples for which MPBNGC was able to terminate successfully, MPBNGC 
was the fastest algorithm again. Considering the competitors, for the nonlinearly con¬ 
strained problems with T = 1 the reduced algorithm is 13.3 seconds resp. 11.3 seconds 
faster than SolvOpt, while for the nonlinearly constrained problems with T = \y\^ 
SolvOpt is 7.1 seconds resp. 5.4 seconds faster than the reduced algorithm. 

Taking a closer look at pi yields the observation that at least 85% of the time 
is consumed by solving the QP (in the linearly constrained case) resp. at least 80% 
of the time is consumed by solving the QCQP (in the nonlinearly constrained case), 
which implies that the difference in the percentage between the QP and the QCQP is 
small in particular (an investigation of the behavior of the solving time ti for higher 
dimensional problems can be found in Fendl & Schichl 0). 

Therefore, we will concentrate in Subsection l4.3l Subsection l4.4l and Subsection l4.5l 
on the comparison of qualitative aspects between the second order bundle algorithm, 
MPBNGC and SolvOpt (like, e.g., the cost c of the evaluations), where before making 
these detailed comparisons, we give a short overview of them as a reason of clarity of 
the presentation: In both cases T = 1 (solid line) and T = \y\^ (dashed line), where we 
use the two different line types for a better distinction in the following, we tested 128 
linearly constrained examples with a fixed box, 117 linearly constrained examples with 
a variable box and 201 nonlinearly constrained examples, which yields the following two 
summary tables consisting of the number of examples for which the second order bundle 
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algorithm (BNLC resp. the reduced algorithm) is better than MPBNGC resp. SolvOpt 
(and vice versa) with respect to the cost c of the evaluations 


(Color code: Light grey) 

MPBNGC 

no termi- significantly better a bit 

nation better better 

nearly 

equal 

BNLC/Red Alg 

a bit better significantly 

better better 

Linearly constrained (fixed box) 
Linearly constrained (variable box) 
Nonlinearly constrained 

T = 1 

0 2 5 12 106 2 0 1 

0 0 0 1 116 0 0 0 

32 6 28 89 31 10 2 3 

Linearly constrained (fixed box) 
Linearly constrained (variable box) 
Nonlinearly constrained 

T = \y\2 

0 2 5 30 91 0 0 0 

0 0 0 5 112 0 0 0 

43 4 28 59 30 15 14 8 

(Color code: Black) 

SolvOpt 

no termi- significantly better a bit 

nation better better 

nearly 

equal 

BNLC/Red Alg 

a bit better significantly 

better better 

Linearly constrained (fixed box) 
Linearly constrained (variable box) 
Nonlinearly constrained 

T = 1 

0 1 3 0 61 25 13 25 

0 0 0 0 48 37 24 8 

0 0 14 20 21 76 20 50 

Linearly constrained (fixed box) 
Linearly constrained (variable box) 
Nonlinearly constrained 

T = \yh 

0 1 2 1 32 34 49 9 

0 0 0 5 41 32 19 20 

0 2 24 26 31 61 45 12 


that are visualized in Figures fT^ 1161 and 1 171 


120 
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Fig. 15: Linearly constrained with fixed box (summary) 
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Fig. 16: Linearly constrained with variable box (summary) 
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Fig. 17; Nonlinearly constrained (summary) 

and that let us draw the following conclusions: 

The performance differences between BNLC and MPBNGC can be neglected for the 
largest part of the linearly constrained examples (with small advantages for MPBNGC 
in about ten percent of these examples). For the nonlinearly constrained examples the 
reduced algorithm is superior to MPBNGC in one quarter of the examples, for forty 
percent of the examples one of these two solvers has small advantages over the other (in 
most cases MPBNGC is the slightly more successful one), the performance differences 
between the two algorithms considered can be completely neglected for fifteen percent 
of the examples, and for further fifteen percent of the examples MPBNGC beats the 
reduced algorithm clearly. 

For the linearly constrained examples BNLC is superior to SolvOpt in one third 
of the examples, for one quarter of the examples one of these two solvers has small 
advantages over the other (in nearly all cases BNLC is the slightly more successful one), 
the performance differences between the two algorithms considered can be completely 
neglected for forty percent of the examples, and in only one percent of the examples 
SolvOpt beats the reduced algorithm clearly. For the nonlinearly constrained examples 
the reduced algorithm is superior to SolvOpt in one third of the examples, for 45 
percent of the examples one of these two solvers has small advantages over the other 
(the reduced algorithm is often the slightly more successful one), the performance 
differences between the considered two algorithms can be completely neglected for ten 
percent of the examples, and in the remaining ten percent of the examples SolvOpt 
beats the reduced algorithm clearly. 

In contrast to the linearly constrained case, in which all three solvers terminated 
successfully for all examples, only the reduced algorithm and SolvOpt were able to 
attain this goal in the nonlinearly constrained case, too. 


4.3 Linearly constrained case (fixed box) 

We took 310 examples from real CSPs that occur in GloptLab. We observe that for 79 
examples the starting point is feasible for the CSP and for 103 examples the evaluation 
of the certificate at the starting point identifies the box as infeasible and hence there 
remain 128 test problems. 

BNLC vs. MPBNGC In the case T = 1 we conclude from Figure fTHl that BNLC is 
significantly better in 1 example and a bit better in 2 examples in comparison with 
MPBNGC, while MPBNGC is significantly better in 2 examples, better in 5 examples 
and a bit better in 12 examples in comparison with BNLC. In the 106 remaining 
examples the costs of BNLC and MPBNGC are practically the same. 

In the case T = lyl^ it follows from Figure IT^ that MPBNGC is significantly better 
in 2 examples, better in 5 examples and a bit better in 30 examples in comparison with 
BNLC. In the 91 remaining examples the costs of BNLC and MPBNGC are practically 
the same. 

BNLC vs. SolvOpt In the case T = 1 we conclude from Figure [20] that BNLC is signif¬ 
icantly better in 25 examples, better in 13 examples and a bit better in 25 examples in 
comparison with SolvOpt, while SolvOpt is significantly better in 1 example and better 
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in 3 examples in comparison with BNLC. In the 61 remaining examples the costs of 
BNLC and SolvOpt are practically the same. 

In the case T = \y\^ it follows from Figure [5T] that BNLC is significantly better in 
9 examples, better in 49 examples and a bit better in 34 examples in comparison with 
SolvOpt, while SolvOpt is significantly better in 1 example, better in 2 examples and 
a bit better in 1 example in comparison with BNLC. In the 32 remaining examples the 
costs of BNLC and SolvOpt are practically the same. 


4.4 Linearly constrained case (variable box) 

We observe that for 80 examples the starting point is feasible for the CSP and for 113 
examples the evaluation of the certificate at the starting point identifies the boxes as 
infeasible and hence there remain 117 test problems of the 310 original examples from 
GloptLab. 

BNLC vs. MPBNGC In the case T = 1 we conclude from Figure \n\ that MPBNGC 
is a bit better in 1 example in comparison with BNLG. In the 116 remaining examples 
the costs of BNLG and MPBNGG are practically the same. 

In the case T = \y\.^ it follows from Figure that MPBNGG is a bit better in 5 
examples in comparison with BNLG. In the 112 remaining examples the costs of BNLG 
and MPBNGC are practically the same. 

BNLC vs. SolvOpt In the case T = 1 we conclude from Figure [24] that BNLC is 
significantly better in 8 examples, better in 24 examples and a bit better in 37 examples 
in comparison with SolvOpt. In the 48 remaining examples the costs of BNLC and 
SolvOpt are practically the same. 

In the case T = \y\.^ it follows from Figure [2^ that BNLC is significantly better 
in 20 examples, better in 19 examples and a bit better in 32 examples in comparison 
with SolvOpt, while SolvOpt is a bit better in 5 examples (21, 101, 102, 128, 189) in 
comparison with BNLC. In the 41 remaining examples the costs of BNLC and SolvOpt 
are practically the same. 


4.5 Nonlinearly constrained case 

Since we were not able to find a starting point, i.e. an infeasible sub-box, for 109 
examples, we exclude them from the following tests for which there remain 201 examples 
of the 310 original examples from GloptLab. 

Reduced algorithm vs. MPBNCC In the case T = 1 MPBNGC does not satisfy any 
of its termination criterions for 32 examples within the given number of iterations 
and function evaluations. For the remaining 169 examples we conclude from Figure !^ 
that the reduced algorithm is significantly better in 3 examples, better in 2 examples 
and a bit better in 10 examples in comparison with MPBNGC, while MPBNGC is 
significantly better in 6 examples, better in 28 examples and a bit better in 89 examples 
in comparison with the reduced algorithm, and in 31 examples the costs of the reduced 
algorithm and MPBNGC are practically the same. 
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In the case T = \y\^ MPBNGC does not satisfy any of its termination criterions 
for 43 examples within the given number of iterations and function evaluations. For 
the remaining 158 examples it follows from Figure [57] that the reduced algorithm is 
significantly better in 8 examples, better in 14 examples and a bit better in 15 examples 
in comparison with MPBNGC, while MPBNGC is significantly better in 4 examples, 
better in 28 examples and a bit better in 59 examples in comparison with the reduced 
algorithm, and in 30 examples the costs of the reduced algorithm and MPBNGC are 
practically the same. 


Reduced algorithm vs. SolvOpt In the case T = 1 we conclude from Figure [5H] that the 
reduced algorithm is significantly better in 50 examples, better in 20 examples and a 
bit better in 76 examples in comparison with SolvOpt, while SolvOpt is better in 14 
examples and a bit better in 20 examples in comparison with the reduced algorithm. 
In the 21 remaining examples the costs of the reduced algorithm and SolvOpt are 
practically the same. 

In the case T = \y\.^ it follows from Figure[5^that the reduced algorithm is signifi¬ 
cantly better in 12 examples, better in 45 examples and a bit better in 61 examples in 
comparison with SolvOpt, while SolvOpt is significantly better in 2 examples, better 
in 24 examples and a bit better in 26 examples in comparison with the reduced algo¬ 
rithm. In the 31 remaining examples the costs of the reduced algorithm and SolvOpt 
are practically the same. 


5 Conclusion 

In this paper we presented a nonsmooth function that can be used as a certificate 
of infeasibility that allows the identification of exclusion boxes during the solution 
process of a GSP by techniques from nonsmooth optimization: While we can find an 
exclusion box by solving a linearly constrained nonsmooth optimization problem, the 
enlargement of an exclusion box can be achieved by solving a nonlinearly constrained 
nonsmooth optimization problem. Furthermore, we discussed important properties of 
the certificate as the reduction of scalability and we suggested a method to obtain a 
good starting point for the nonsmooth optimization problems. 
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Fig. 18; Linearly constrained (fixed box) — rp(c) for BNLC & MPBNGC (T = 1) 
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Fig. 19: Linearly constrained (fixed box) — rp(c) for BNLC Sz MPBNGC (T = \y\ 2 ) 
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Fig. 20: Linearly constrained (fixed box) — rp(c) for BNLC & SolvOpt (T = 1) 
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Fig. 21: Linearly constrained (fixed box) 


rp(c) for BNLC SolvOpt (T 
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Fig. 22: Linearly constrained (variable box) — rp(c) for BNLC & MPBNGC (T = 1) 
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Fig. 23: Linearly constrained (variable box) — rp(c) for BNLC & MPBNC (T = I 2 /I 2 ) 
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Fig. 24: Linearly constrained (variable box) — rp(c) for BNLC Sz SolvOpt (T = 1) 
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Fig. 26: Nonlinearly constrained — rp(c) for Red Alg & MPBNGC (T = 1) 
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Fig. 27: Nonlinearly constrained — rp(c) for Red Alg & MPBNGC (T = 
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Fig. 28: Nonlinearly constrained — rp(c) for Red Alg & SolvOpt (T = 1) 
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Fig. 29: Nonlinearly constrained — rp(c) for Red Alg & SolvOpt (T = \y\ 2 ) 


















